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Using ab initio lattice methods, we calculate the finite temperature thermodynamics of homo¬ 
geneous two-dimensional spin-1/2 fermions with attractive short-range interactions. We present 
results for the density, pressure, compressibility, and quantum anomaly (i.e. Tan’s contact) for a 
wide range of temperatures and coupling strengths, focusing on the unpolarized case. Within our 
statistical and systematic uncertainties, our prediction for the density equation of state differs quan¬ 
titatively from the prediction by Luttinger-Ward theory in the strongly coupled region of parameter 
space, but otherwise agrees well with it. We also compare our calculations with the second- and 
third-order virial expansion, with which they are in excellent agreement in the low-fugacity regime. 
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Introduction.- Quantum mechanics in two spatial di¬ 
mensions (2D) is a fascinating playground to understand 
fundamental physics in a wide variety of situations. It 
represents a necessary (though often odd!) crossover be¬ 
tween the integrable systems that live in ID and their 
much more challenging 3D counterparts. Interest in 2D 
ranges from basic theory and experiment to marketable 
technological applications. Among the most salient fea¬ 
tures and systems in 2D we have: the peculiar behavior 
of Kosterlitz-Thouless phase transitions [I] ; the possibil¬ 
ity of understanding quark confinement analytically in 
certain models [2]; anomalous scale invariance in non- 
relativistic systems !3]; the challenge of high-temperature 
superconductors and, of course, graphene [5]. 

In recent years, experiments with ultracold atoms HE] 
have made clear progress towards a clean and system¬ 
atic study of fermionic atom clouds in constrained or 
modulated optical traps, which closely resemble a 2D 
situation [M3] (see also Ref. m for a recent discus¬ 
sion approaching that limit). These systems have been 
of singular interest due to their malleability: by now 
it is well known that the inter-atomic interaction can 
be tuned essentially at will, low-temperature degenerate 
regimes can be reached, spin- and mass-asymmetric sys¬ 
tems can be studied, and so on. As advances are thus 
made towards understanding the thermodynamics, struc¬ 
ture and phases on the experimental side, theorists are 
once again faced with the challenge of strongly coupled 
regimes, which defy analytic approaches. 

Previous theoretical studies have characterized the 
crossover between Bose-Einstein condensation (BEC) 
and BCS pairing in 2D via mean-field analyses [151 - 
E23- The first determination of the ground-state equa¬ 
tion of state, however, appeared relatively recently in 
Ref. [T5], which used the diffusion Monte Carlo method. 
Even more recently, Ref. [1DJ performed an auxiliary-field 
quantum Monte Carlo study of the ground state where 
the pressure, contact, and condensate fraction were de¬ 
termined. At finite temperature, the equation of state 
was first computed in the virial expansion in Ref. m, 


and in the Luttinger-Ward approach in Ref. [5Tj. Pair 
correlations were investigated using the virial expansion 
at about the same time in Ref. [22] and in Ref. [23]. Col¬ 
lective modes were studied in Refs. [24lt26j . and the shear 
viscosity and spin diffusion in Ref. [ 13 - 

In this work we offer a few essential pieces of the ther¬ 
modynamic puzzle of attractively interacting fermions in 
2D. We implement a lattice Monte Carlo method to cal¬ 
culate the thermal and short-range correlation properties 
of the system along the 2D analogue of the BCS-BEC 
crossover. Specifically, we determine the density equa¬ 
tion of state, from which we extract the pressure and 
compressibility; all of these are experimentally measur¬ 
able. To our knowledge, the present is the first fully 
non-perturbative, ab initio calculation of the equation of 
state of this system at finite temperature. In addition, 
we present here the first calculation of Tan’s contact pa¬ 
rameter [28H30] at finite temperature in 2D that is free 
of uncontrolled approximations. 

Hamiltonian and many-body method.- The dynamics 
of dilute spin-1/2 non-relativistic fermions with short- 
range interactions is governed by the following Hamilto¬ 
nian: 


V’sM ^ S’V X )™4-( X ) 

( 1 ) 

where ?/*], ^ s are the creation and annihilation operators 
in coordinate space for spin s particles, and h s = ’4>14> S 
are the corresponding densities. 



In a recent work [31) . we performed a similar study to 
the one presented here, but restricted to ID. Here, we use 
the same methods (which are also closely related to those 
of Refs. pT2TPm used in 3D) but applied in 2D. We placed 
the spin-1/2 Fermi system in a Euclidean spacetime lat¬ 
tice of extent N% x N T with periodic boundary conditions 
in the spatial directions and anti-periodic in the time di¬ 
rection. A Trotter-Suzuki factorization, followed by a 
Hubbard-Stratonovich transformation [55] , allowed us to 
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write the grand-canonical partition function as a path 
integral over an auxiliary field. That path integral was 
evaluated using Metropolis-based Monte Carlo methods 
(see, e.g., Refs. 37]). Throughout this work, we use 
units such that h = m = k-Q = 1, where m is the mass of 
the fermions. The physical spatial extent of the lattice is 
Lx L, where L = N x t and we take £ = 1 to set the length 
and momentum scales. The extent of the temporal lat¬ 
tice is set by the inverse temperature /3 = 1/T = tN t . 
The time step r = 0.05 (in lattice units) was chosen 
to balance temporal discretization effects with computa¬ 
tional efficiency; in any case, those discretization effects 
are smaller than our statistical effects. 

As in our previous study, the physical input parameters 
are the inverse temperature /3, the chemical potential g = 
/T* = g, , and the (bare, attractive) coupling strength g > 
0, where the latter determines the binding energy e B of 
the two-body problem (see below). From these, we form 
two dimensionless quantities: the fugacity 2 = exp(/3/i) 
and the dimensionless coupling /3e B . Using z and j3e B 
as parameters facilitates comparisons with experiments, 
as well as with other theoretical approaches such as the 
Luttinger-Ward calculations mentioned above. 

The connection between the bare coupling g and the 
2D scattering length a 2D via the two-body problem re¬ 
quires us to set the scales £ and L. To avoid ambigui¬ 
ties, we computed the binding energy e B of the two-body 
problem on each of the lattices studied, which allows us 
to characterize the coupling in terms of the dimension¬ 
less parameter (3e b above. This 2D case, in contrast to 
its ID and 3D analogues, is peculiar due to its anomalous 
scale invariance. Indeed, as may be inferred from Eq. ([!]), 
the dimensions of if a are of inverse length (in the natu¬ 
ral units mentioned above), such that g is dimensionless 
and therefore the system is classically (i.e. before ac¬ 
counting for quantum fluctuations) scale-invariant. The 
appearance of a two-body bound state with binding en- 
ergy e B , upon turning on the interaction, is an example 
of a quantum anomaly: the classical scale invariance of 
H is broken by quantum effects. Notably, massless quan¬ 
tum chromodynamics in 3D, also a scale-free theory at 
the classical level, displays a similar feature. 

Lattice calculations of the kind we use are exact, up to 
statistical and systematic uncertainties. To address the 
former, we took 1000 de-correlated samples for each data 
point in the plots shown below, which yields a statistical 
uncertainty of order 3%. The smoothness of our results 
show that those effects are indeed small. 

To address the systematic effects, one must approach 
the continuum limit. Two-dimensional problems are 
computationally inexpensive relative to full 3D problems, 
but they are still challenging. In this first study, we re¬ 
stricted ourselves to N x = 11,15,19. The continuum 
limit is achieved by lowering the density while still re¬ 
maining in the many-particle, thermodynamic regime. In 
turn, this is accomplished by increasing the lattice pa- 



Figure 1. (color online) Density equation of state, in units of 
the non-interacting density n 0 of spin-1/2 fermions in 2D, for 
coupling strengths /3e B = 0.1, 0.5, 1, 2, 3 (from bottom to 
top), as a function of fig. The error bars reflect the statistical 
uncertainty. The solid colored lines show the Luttinger-Ward 
result of Ref. m- The long- and short-dashed lines show the 
second- and third-order virial expansion results, respectively. 

rameter /3, ensuring that the thermal wavelength X T = 
\/2i t(3 satisfies £ < A T < L; at fixed z, this reduces 
the density. In this work, we used X T ~ 8.0, which is 
well within the regime specified above. We then studied 
whether our results collapse to a universal curve when /3 
and g are varied while (3e B is held fixed. Lattice sizes 
larger than N x = 19 are computationally more expen¬ 
sive but not impractical; however, we chose to fix that 
size and cover a wider region of parameter space instead. 
Because our study proceeded at constant j3e B , increasing 
/3 implies reducing g , which results in smaller uncertain¬ 
ties associated with the temporal lattice spacing r in the 
Trotter-Suzuki decomposition; these are expected to be 
of order 1 — 2% (see e.g. Ref. [33]). 

Analysis and Results .- To characterize the thermo¬ 
dynamics of a strongly interacting system, as is our ob¬ 
jective here, a simple yet extremely effective route is to 
first calculate the total particle number density N/L 2 
(where N = N^. + and N s is the particle number 
for spin s =f, 4-) as a function of the thermodynamic 
variables (temperature, chemical potential, and interac¬ 
tion strength, as mentioned above). The density has the 
added benefit over other quantities that its statistical- 
noise effects in a Monte Carlo calculation are relatively 
small. From n, one may determine the pressure P by 
performing a numerical integration along fig, and the 
compressibility n by differentiation. 

In Fig. [l] we show the density n as a function of the 
dimensionless parameters 2 and /3e B , defined above. The 
non-interacting result used as a scale in that figure is 
n 0 X t = 4 ^( 2 ), where Ii(z) = zdl 0 (z)/dz, and 

/»oo 

I 0 (z) = / dx xki(l + ze~ x ). ( 2 ) 

Jo 
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Our calculations, as shown in Fig. [T] display a clear 
deviation from the Luttinger-Ward approach of Ref. |2T] 
for —2 < j3fi < 2. This is not unexpected, as that regime 
is challenging for most many-body approaches: it is the 
regime where quantum fluctuations are the strongest (at 
fixed /3e B ). Away from that region, however, the agree¬ 
ment with Ref. 12H satisfactory. Moreover, both cal¬ 
culations seem to heal together to the virial expansion 
result at large and negative py. The starting point for 
this expansion is a Taylor expansion of the grand ther¬ 
modynamic potential Q in powers of the fugacity z, 


— /3fl — Q 1 (z + 6 2 z 2 + 6 3 z 3 + ...) , (3) 

where Q ± = 2L 2 /\\ is the 2D single-particle partition 
function and b n are the virial coefficients, which in gen¬ 
eral will be functions of the coupling pe B . Thus, the virial 
expansion for the density reads nX^/2 = z + 2b 2 z 2 + 
3 b 3 z 3 + • • •, where the factor of 1/2 on the left-hand 
side comes from the number of fermion species. For the 
system studied here, the second-order coefficient b 2 is 
known analytically from the exact solution of the two- 
body problem (see e.g. Refs. [221I2T| ): 


b, = -- 




dy 2e 


y 7r 2 + 4 hr y 


(4) 


A calculation of b 3 can be found in Ref. [22] , 

To calculate the pressure P, we integrate as follows 


P\^ = 2t r 



(5) 


where we have put everything in dimensionless form us¬ 
ing the thermal wavelength scale. In Fig. [2] we show P, 



Py 



Figure 3. (color online) Compressibility, in units of the non¬ 
interacting compressibility k 0 , of spin-1/2 fermions in 2D for 
coupling strengths /3e b = 0.1, 0.5, 1, 2, 3 (from top to bot¬ 
tom), as a function of py. The error bars reflect the dif¬ 
ference between a smooth interpolation and the raw Monte 
Carlo data. 


as obtained from the above formula, in units of the pres¬ 
sure of the non-interacting system P 0 A^ = 8ttI q (z). The 
virial expansion of Eq. ([3]) is used in this integration to 
complete the approach to the z —> 0 limit. 

On the other hand, by taking a derivative of n one 
obtains the isothermal compressibility, 


K 


• dn 


d(Py) 




d{n\\) 


27t (n\\) 2 d(py) 


( 6 ) 


We report this quantity in Fig. [3j in units of its non¬ 
interacting counterpart n 0 , where (in dimensionless form) 
k 0 Ap 4 = (2/7r)(n 0 Ay) _2 / 2 (z), and / 2 (z) = zdl 1 (z)/dz. 

To calculate the contact, we use the grand-canonical 
definition 1221 




d(psi) 


P <91n(a 2D /A T ) 


(7) 


T,n 


From here, it is easy to see that the virial expansion for 
C takes the form 


PC = 2nQ 1 (c 2 z 2 + c 3 z 3 + ...), (8) 

where c n = —db n /d ln(a 2D /A T ). Using Eq. [4j it is 
straightforward to obtain the exact continuum-limit an¬ 
swer: 


C 2 


2pe B e f5s K 


1 + 2 



y e -^ B (2/ 2 + l) 

d V— 2 - 75 2 

7r 2 + 4 In y 


(9) 


Figure 2. (color online) Pressure, in units of the non¬ 
interacting pressure P 0 , of spin-1/2 fermions in 2D for cou¬ 
pling strengths /?e B = 0.1, 0.5, 1, 2, 3 (from bottom to top), 
as a function of Py. The error bars reflect the statistical un¬ 
certainty. The long- and short-dashed lines show the second- 
and third-order virial expansion result, respectively. 


This result was evidently used in Ref. [2^, but the for¬ 
mula itself is not shown explicitly anywhere, to the best 
of our knowledge. 

In our Monte Carlo calculations, we determine the con¬ 
tact by calculating the expectation value of the interac¬ 
tion energy V. Using the definition of Eq. [ 7 J along with 
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— /3Q = In Z , where Z is the grand-canonical partition 
function, we obtain 

-2tt 5 In Z _ n - ding 
P 31n(a 2D /A T ) “ ' 7 dln(a 2D /A T ) ’ 1 j 

where we have used the fact that V is a contact inter¬ 
action as in Eq. [I] The remaining factor on the right 
is determined by solving the 2-body problem. To turn 
this into an intensive, dimensionless quantity, we present 
results in the form C/{Nk 2 F ). In Fig. [ 4 ] we show our 
lattice Monte Carlo results for this quantity as a func¬ 
tion of T/Tp and the coupling f3e B , as well as selected 
lines of constant ln(fc F a 2 D)- Corresponding to the lat¬ 
ter, we have included ground-state quantum Monte Carlo 
(QMC) results [T8] and experimental results at finite tem¬ 
perature of T/T F = 0.27 [IT] . The experimental results 
largely agree with our calculations, however, the experi¬ 
mental errors and the maximum coupling calculated limit 
us from drawing any strong conclusions. On the other 
hand, we note that the contact is largely flat at constant 
ln(fcj,a 2 D)j which agrees qualitatively with the Luttinger- 
Ward approach of Ref. m ■ We regard all this as evi¬ 
dence that the ground-state QMC results (see also next 
section) are in very good agreement with our finite-T 
calculations, as the latter seem to approach the T = 0 
results in that limit. 

Our full set of contact calculations as a function of /3/jl 
can be found in the Supplemental Material. We have 
also verified there the agreement with the continuum- 
limit second-order virial expansion result (using Eq. [ 9 ]) 
in the high-energy regime. Removing lattice effects in 
this regime is most demanding, which suggests that such 
systematic effects are well under control. 

Additional Tests and Crosschecks.- As a test, 
we used the same lattice approach, in combination 
with the projection Monte Carlo method, to calcu¬ 
late the ground-state energies of the system as a func¬ 
tion of ln(fc F a 2D ), which were first determined non- 
perturbatively in Ref. [18] , Our results, along with those 


Table I. Ground-state energy in units of its non-interacting 
counterpart, as a function of the coupling ln(fc F a 2D ). Here, 
a 2D = a 2 D e 7 /2, where 7 ~ 0.577 is Euler’s constant, which is 
the convention followed in Ref. [18l . 


ln(fc F fi 2D ) 

Ags/^fg 1 Ref- I18j 

This work 

Difference 

- 2.00 

-137.761(7) 

-139.4(5) 

1 % 

-1.50 

-50.593(4) 

-51.0(8) 

1 % 

- 1.00 

-18.532(4) 

— 19.0(5) 

3% 

-0.50 

-6.714(4) 

- 6 . 8 ( 2 ) 

1 % 

0.00 

-2.318(2) 

-2.40(5) 

3% 

0.25 

-1.283(12) 

— 1.35(3) 

5% 

0.50 

—0.638(10) 

-0.70(2) 

10 % 

1.44 

0.349(6) 

0.38(1) 

9% 

3.34 

0.706(2) 

0.698(1) 

1 % 


4 


Nk F 2 

2 


1 


T l T F 

Figure 4. (color online) Contact, in units of l/(Nk%), for 
spin-1/2 fermions in 2D for coupling strengths /3e b = 0.1, 
0.5, 1, 2, 3 (from bottom to top), as a function of T/Tf. 
The error bars reflect the statistical uncertainty. Solid lines 
of gray across colored contact results indicate lines of con¬ 
stant ln(fc F a 2 D) = 1.60,1.21,0.85,0.53,0.24 from bottom to 
top. Corresponding experimental data at T/Tf = 0.27 from 
Ref. Em given by solid circles (displaced about light gray 
line for visibility of error bars); corresponding Luttinger-Ward 
results at T/Tf = 0.27 from Ref. ( 2 T] are shown as black 
stars; and corresponding QMC calculations from Ref. m at 
T/Tf = 0 given by solid squares and dotted gray lines. 



of Ref. (IE] , are shown in Table |T] Note that Ref. [TS] 
uses ln(fcpd 2D ) as the coupling, where a 2D = a 2 D e 7 /2 
and 7 ~ 0.577 is Euler’s constant. 

Summary and conclusions.- Using lattice Monte Carlo 
methods, we have calculated the finite temperature ther¬ 
modynamics of homogeneous two-dimensional spin- 1/2 
fermions with attractive short-range interactions. We 
have presented results for the density, pressure, com¬ 
pressibility, and Tan’s contact for a wide range of temper¬ 
atures and coupling strengths. Within our statistical and 
systematic uncertainties, our prediction for the density 
equation of state differs from the prediction by Luttinger- 
Ward theory in a substantial region of parameter space. 
The general agreement is nonetheless exceptional. We 
have also compared our calculations of the density and 
pressure with the second- and third-order virial expan¬ 
sion, with which they agree remarkably well in the low 
fugacity regime. Moreover, the agreement seems stronger 
with our results than with the Luttinger-Ward approach. 
Finally, we have presented a comparison of our calcula¬ 
tion of the contact with previous ground-state calcula¬ 
tions and finite-temperature experimental data. A more 
complete representation of our data for the contact, in¬ 
cluding a comparison with the second-order virial expan¬ 
sion and an alternative temperature scale, appear in the 
Supplemental Material. 

Our results for the density, pressure, and compressibil¬ 
ity can also be compared with experiments. One of the 
motivations for the latter is that attractively interact- 
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ing fermions in 2D are expected to undergo a Kosterlitz- 
Thouless transition into a superfluid phase at low enough 
temperatures. We do not see, in the quantities studied 
here, any particular signature of the transition. We defer 
further calculations in that direction to future work. 
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SUPPLEMENTAL ONLINE MATERIAL 


We display our data for the temperature and contact 
in scales different from those in the main text. We com¬ 
pare the contact explicitly with the second-order virial 
expansion. 



Pn 


Figure 6. (color online) Contact of spin-1/2 fermions in 2D, 
for coupling strengths /3e b = 0.1, 0.5, 1, 2, 3 (from bottom to 
top), as a function of /3/u. 


TEMPERATURE SCALE 


Although fin is a useful dimensionless parameter to 
characterize the temperature, especially in the grand- 
canonical ensemble, it is often even more useful to present 
the temperature in a scale set by the density, namely 
the Fermi energy e F = k 2 F / 2, where k F = y/2nn is the 
Fermi momentum. In Fig. [5] we show T/e f as a func¬ 
tion of ln(/c F a 2D ) = ln(2e F /e B )/2, where a 2D = 

The statistical uncertainties are not shown, but from the 
smoothness of the resulting curves it can be inferred that 
they are very small. Their size results from the tiny un¬ 
certainty in the density n (see density plot in main text), 
which (barely) impacts both axes in Fig. [5]. 



Figure 5. (color online) Temperature scale T, in units of the 
non-interacting Fermi energy e F , of spin-1/2 fermions in 2D, 
for coupling strengths /3e b = 0.1, 0.5, 1, 2, 3 (from top to 
bottom), as a function of the more common coupling strength 
parameter ln(fc F a 2D ). The dashed line is a fit to the /3/x = 0 
points for each coupling; the fit form is A(\n(k F a 2 Y)) — 0.5) s , 
with A = 1.170 and B = 0.422. The value at T = 0, namely 
ln(fc F a 2D ) — 0.5, was estimated from Ref. [Tj. 


TAN’S CONTACT 

In the main text we showed Tan’s contact C as a func¬ 
tion of T/e F . However, in our calculations e F is com¬ 
puted as an output after the average density is deter¬ 
mined. Therefore, it is more natural to report C as a 
function of /3/x, which we do in Fig. [6] The statistical 
uncertainties in our Monte Carlo data are smaller than 
the size of the symbols in that figure. 

An additional benefit of the above representation is 
that it allows for a more direct comparison with the virial 
expansion at large and negative /3/z. Our calculations 
agree with that expansion to second order (leading order 
for the contact), as can be appreciated in Fig. [7] The 
statistical uncertainties can be seen in this case. 

Finally, to further complement the contact plots of the 
main text, we provide in Fig. [8]the full temperature range 
we explored in our calculations. 
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Figure 7. (color online) Contact of spin-1/2 fermions in 2D, 
for coupling strengths /3e B = 0.1, 0.5, 1, 2, 3 (from bottom 
to top), as a function of /3/j, in the low-fugacity regime. The 
dashed lines show the second-order virial expansion result. 
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Figure 8. (color online) Contact of spin-1/2 fermions in 2D 
as a function of /3/r. The solid colored lines show lines of 
constant /3s B = 0.1, 0.5, 1, 2, 3 (from bottom to top). Solid 
lines of gray across colored contact results indicate lines of 
constant ln(fc F a 2D )=l-25, 1.00, 0.75, 0.5, 0.25, 0, -0.25, -0.5 
from bottom to top. 















